TMB of tumors across unpaired longitudinal samples in the PBTA Cohort

Tumor evolution project

Data used

In this notebook, we are using the pbta.tsv file geenrated by the sample-distribution-analysis module.

Set up

suppressPackageStartupMessages({
  library(tidyverse)
  library(rstatix)
  library(ggpubr)
})

Directories and File Inputs/Outputs

# Detect the ".git" folder -- this will be in the project root directory
# Use this as the root directory to ensure proper sourcing of functions 
# no matter where this is called from
root_dir <- rprojroot::find_root(rprojroot::has_dir(".git"))
analysis_dir <- file.path(root_dir, "analyses", "tmb-vaf-longitudinal")
input_dir <- file.path(analysis_dir, "input")
files_dir <- file.path(root_dir, "analyses", "sample-distribution-analysis", "results")

# Input files
pbta_file <- file.path(files_dir, "pbta.tsv") # file from add-sample-distribution module
tmb_file <- file.path(input_dir, "snv-mutation-tmb-coding.tsv")
palette_file <- file.path(root_dir, "figures", "palettes", "tumor_descriptor_color_palette.tsv")

# File path to plot directory
plots_dir <-
  file.path(analysis_dir, "plots")
if (!dir.exists(plots_dir)) {
  dir.create(plots_dir)
}

source(paste0(root_dir, "/figures/scripts/theme.R"))

Read in data and process

# Make this reproducible
set.seed(2023)

# Read in tmb file and remove hypermutant samples
tmb <- readr::read_tsv(tmb_file, guess_max = 100000, show_col_types = FALSE) %>%
  filter(!tmb >= 10) %>% 
  dplyr::rename(Kids_First_Biospecimen_ID = Tumor_Sample_Barcode) # change name of the biospecimen to match the one from the histologies files

# Read in pbta file and create df_plot
df_plot <- readr::read_tsv(pbta_file, guess_max = 100000, show_col_types = FALSE) %>%
  filter(!(experimental_strategy == "RNA-Seq"),
         !(tumor_descriptor == "Unavailable")) %>% # these are 2 samples
  left_join(tmb, by = c("Kids_First_Biospecimen_ID", "experimental_strategy")) %>%
  select(Kids_First_Participant_ID, Kids_First_Biospecimen_ID, tumor_descriptor,
         broad_histology, match_id_DNA, experimental_strategy, 
         cancer_group, short_histology, tmb, mutation_count, DNA) %>%
  dplyr:::mutate(short_histology_sum = ifelse(short_histology == "HGAT", "HGG",
                                              ifelse(short_histology == "LGAT", "LGG", 
                                                     ifelse(short_histology == "Ependymoma", "Ependymoma", 
                                                            ifelse(short_histology == "Medulloblastoma", "Medulloblastoma", "Other"))),
                 log10_tmb = abs(log10(tmb)),
                 tumor_descriptor = factor(tumor_descriptor),
                 tumor_descriptor = fct_relevel(tumor_descriptor, c("Diagnosis", "Progressive", "Recurrence", "Deceased", "Second Malignancy", "Unavailable"))) %>% 
  distinct(match_id_DNA, .keep_all = TRUE) %>% # keep only one bs_sample per genomic assay (WGS/WXS)
  filter(!is.na(tmb))

# How many samples per Timepoint?
print(df_plot %>% count(tumor_descriptor))

Error: unexpected symbol in:
"# How many samples per Timepoint?
print"

How many patient samples in the unpaired longitudinal data?

There are 1924 samples in total.

All unpaired samples

We are interested in how TMB changes over the time in all PBTA samples (unpaired samples).

# Read color palette
palette_df <- readr::read_tsv(palette_file, guess_max = 100000, show_col_types = FALSE) %>% 
  mutate(tumor_descriptor = color_names)

# Define and order palette
palette <- palette_df$hex_codes
names(palette) <- palette_df$tumor_descriptor

# Define label for plots
Timepoint <- df_plot$tumor_descriptor

# Define ylim
ylim <- max(df_plot$tmb)

print(df_plot %>% 
        ggplot2::ggplot() + 
        ggplot2::aes(x = tumor_descriptor, 
                     y = tmb,
                     color = Timepoint) + 
        ggplot2::geom_jitter(width = 0.2, alpha = 0.5, size = 1.5) + 
        # light guiding line representing 0 exposure
        ggplot2::geom_hline(yintercept = 0, size = 0.15) + 
        ggplot2::scale_color_manual(values = palette, 
                                    breaks = sort(names(palette))) +
        ggplot2::stat_summary(color = "black", size = 0.3) + 
        ggplot2::labs(x = "Tumor descriptor",
                      y = "TMB") +
        theme_Publication() +
        scale_y_continuous(limits = c(0, ylim)) +
        theme(axis.text.x = element_text(angle = 90)))

# Save the plot
ggsave(filename = "TMB-unpaired-longitudinal-samples.pdf",
       path = plots_dir, 
       width = 10, 
       height = 5, 
       device = "pdf", 
       useDingbats = FALSE)

All unpaired samples across cancer types

We will use the short_histology column to display TMB differences across cancer types.

print(df_plot %>%
        ggplot2::ggplot() + 
        ggplot2::aes(x = tumor_descriptor, 
                     y = tmb,
                     color = Timepoint) + 
        ggplot2::geom_jitter(width = 0.2, alpha = 0.5, size = 1.5) + 
        # light guiding line representing 0 exposure
        ggplot2::geom_hline(yintercept = 0, size = 0.15) + 
        ggplot2::scale_color_manual(values = palette, 
                                    breaks = sort(names(palette))) +
        ggplot2::stat_summary(color = "black", size = 0.3) + 
        ggplot2::facet_wrap(~short_histology, nrow = 8) +
        ggplot2::labs(x = "Tumor descriptor",
                      y = "TMB") +
        theme_Publication() + 
        theme(axis.text.x = element_text(angle = 90)) +
        scale_y_continuous(limits = c(0, ylim)))

# Save the plot
ggsave(filename = "TMB-unpaired-longitudinal-samples-cancer-type.pdf", 
       path = plots_dir, 
       width = 25, 
       height = 20, 
       device = "pdf", 
       useDingbats = FALSE)

We observe that majority of samples are High- and Low-grade gliomas as well as Ependymoma, so we will classify accordingly (versus any other cancer types).

All unpaired samples across cancer types summarized

# Create col with sample size per cancer type
df_plot_n <- df_plot %>%
  dplyr::count(tumor_descriptor) %>%
  dplyr::distinct() %>%
  dplyr::inner_join(
    dplyr::select(palette_df, dplyr::contains("tumor_descriptor"))) %>%
  dplyr::select(tumor_descriptor, n) %>%
  # Create wrapped with (n=X) factor column for cancer groups
  dplyr::mutate(tumor_descriptor_n = glue::glue("{tumor_descriptor} (N={n})")) %>%
  dplyr::inner_join(df_plot) %>%
  unique() %>%
  # reverse order levels to alphabetize when flipping coordinates
  mutate(tumor_descriptor_n = fct_rev(tumor_descriptor_n)) 

Joining with `by = join_by(tumor_descriptor)`Joining with `by = join_by(tumor_descriptor)`

print(df_plot %>% 
        ggplot2::ggplot() + 
        ggplot2::aes(x = tumor_descriptor, 
                     y = tmb,
                     color = Timepoint) + 
        ggplot2::geom_jitter(width = 0.2, alpha = 0.5, size = 1.5) + 
        # light guiding line representing 0 exposure
        ggplot2::geom_hline(yintercept = 0, size = 0.15) + 
        ggplot2::scale_color_manual(values = palette, 
                                    breaks = sort(names(palette))) +
        ggplot2::stat_summary(color = "black", size = 0.3) + 
        ggplot2::facet_wrap(~short_histology_sum, nrow = 1) +
        ggplot2::labs(x = "Tumor descriptor",
                      y = "TMB") +
        theme_Publication() + 
        theme(axis.text.x = element_text(angle = 90)) +
        scale_y_continuous(limits = c(0, ylim)))

# Save the plot
ggsave(filename = "TMB-unpaired-longitudinal-samples-cancer-type-sum.pdf", 
       path = plots_dir, 
       width = 12, 
       height = 6, 
       device = "pdf", 
       useDingbats = FALSE)

Summary statistics for Timepoints and cancer types

Considering the inequalities in the sampling effort for each timepoint, let’s explore whether the sampling size affects the observed patterns, e.g., indication of Progressive and Recurrence having higer TMB compared to Diagnosis for High-grade glioma patients.

# Inspect some random rows of the data by groups
set.seed(2023)
df_plot %>% 
  select(short_histology_sum, tumor_descriptor, tmb, log10_tmb) %>% 
  sample_n_by(short_histology_sum, tumor_descriptor, tmb, size = 1)

# Compute some summary statistics (mean and sd) by groups:
stat.summary <- df_plot %>%
  group_by(tumor_descriptor, short_histology_sum) %>%
  get_summary_stats(tmb, type = "mean_sd") %>% 
  filter(!n <= 2) # remove if total number of values per group is less than 2
stat.summary

# Pairwise comparisons between time points at each group levels
# Paired t-test is used because we have repeated measures by time
stat.test <- df_plot %>%
  group_by(short_histology_sum) %>%
  pairwise_t_test(log10_tmb ~ tumor_descriptor, 
                  paired = FALSE, 
                  p.adjust.method = "bonferroni")
stat.test

# Add statistical test p-values
stat.test <- stat.test %>% add_xy_position(x = "tumor_descriptor")

# Create bxp
print(ggboxplot(df_plot,
                x = "tumor_descriptor",
                y = "log10_tmb",
                color = "tumor_descriptor",
                palette = palette) +
        facet_wrap(~short_histology_sum, nrow = 1) +
        theme_Publication() + 
        theme(axis.text.x = element_text(angle = 90)) +
        stat_pvalue_manual(stat.test, label = "p.adj.signif",
                           step.increase = 0.08, hide.ns = TRUE, tip.length = 0))

# Save the plot
ggsave(filename = "TMB-Bxp-stat-test.pdf", 
       path = plots_dir, 
       width = 12, 
       height = 6, 
       device = "pdf", 
       useDingbats = FALSE)


# Create jitter plot
print(df_plot %>% 
        ggplot2::ggplot() + 
        ggplot2::aes(x = tumor_descriptor, 
                     y = log10_tmb,
                     color = Timepoint) + 
        ggplot2::geom_jitter(width = 0.2, alpha = 0.5, size = 1.5) + 
        # light guiding line representing 0 exposure
        ggplot2::geom_hline(yintercept = 0, size = 0.15) +
        ggplot2::scale_color_manual(values = palette, 
                                    breaks = sort(names(palette))) +
        facet_wrap(~short_histology_sum, nrow = 1) +
        theme_Publication() + 
        theme(axis.text.x = element_text(angle = 90)) +
        stat_pvalue_manual(stat.test, label = "p.adj.signif",
                           step.increase = 0.08, hide.ns = TRUE, tip.length = 0))

# Save the plot
ggsave(filename = "TMB-jitter-stat-test.pdf", 
       path = plots_dir, 
       width = 12, 
       height = 6, 
       device = "pdf", 
       useDingbats = FALSE)

# This code was inspired by the following:
# HOW TO PERFORM MULTIPLE PAIRED T-TESTS IN R
# https://www.datanovia.com/en/blog/how-to-perform-multiple-paired-t-tests-in-r/

sessionInfo()

R version 4.2.3 (2023-03-15)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.4.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] ggthemes_4.2.4  ggpubr_0.6.0    rstatix_0.7.2   lubridate_1.9.2
 [5] forcats_1.0.0   stringr_1.5.0   dplyr_1.1.2     purrr_1.0.1    
 [9] readr_2.1.4     tidyr_1.3.0     tibble_3.2.1    ggplot2_3.4.2  
[13] tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.0  xfun_0.39         bslib_0.5.0       carData_3.0-5    
 [5] colorspace_2.1-0  vctrs_0.6.3       generics_0.1.3    htmltools_0.5.5  
 [9] yaml_2.3.7        utf8_1.2.3        rlang_1.1.1       jquerylib_0.1.4  
[13] pillar_1.9.0      glue_1.6.2        withr_2.5.0       bit64_4.0.5      
[17] lifecycle_1.0.3   munsell_0.5.0     ggsignif_0.6.4    gtable_0.3.3     
[21] ragg_1.2.5        evaluate_0.21     labeling_0.4.2    knitr_1.43       
[25] tzdb_0.4.0        fastmap_1.1.1     parallel_4.2.3    fansi_1.0.4      
[29] highr_0.10        broom_1.0.5       scales_1.2.1      backports_1.4.1  
[33] cachem_1.0.8      vroom_1.6.3       jsonlite_1.8.7    abind_1.4-5      
[37] systemfonts_1.0.4 farver_2.1.1      bit_4.0.5         textshaping_0.3.6
[41] hms_1.1.3         digest_0.6.33     stringi_1.7.12    rprojroot_2.0.3  
[45] cli_3.6.1         tools_4.2.3       magrittr_2.0.3    sass_0.4.7       
[49] crayon_1.5.2      car_3.1-2         pkgconfig_2.0.3   timechange_0.2.0 
[53] rmarkdown_2.23    R6_2.5.1          compiler_4.2.3